﻿using System;
using System.Diagnostics;
using System.Collections.Generic;
using System.Collections;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;
using Emgu.CV;
using System.IO;
using Emgu.CV.Structure;
using Emgu.Util;
using Emgu.CV.UI;


namespace GeometryProcessing
{
    enum ProcedureState
    {
        SUCCESS = 0,
        FAIL = 1
    }

    class MyPointF
    {
        float x; float y;
        public MyPointF(PointF a)
        {
            x = a.X; y = a.Y;
        }
        public MyPointF(float X, float Y) { x = X; y = Y; }
        public static implicit operator PointF(MyPointF p)
        {
            return new PointF(p.x, p.y);
        }
        public static MyPointF operator *(float k, MyPointF a)
        {
            return new MyPointF(a.x * k, a.y * k);
        }
        public static float operator *(MyPointF a, MyPointF b)
        {
            return a.x * b.x + a.y * b.y;
        }
        public static MyPointF operator -(MyPointF a, MyPointF b)
        {
            return new MyPointF(a.x - b.x, a.y - b.y);
        }
    }

    class Gridding
    {
        /*
         * 
         * DEBUG:
         * 
         */
            public bool DEBUG { set; get; }
        /*
         * 
         * DEBUG:
         * 
         */


        public const float EPS = 1e-5F;
        // 划分的数据结构
        // 区域泛型列表
        public ArrayList cellArray { set; get; }
        // 划分的边集合
        public List<LineSegment2DF> lineSegmentArray { set; get; }
        // 划分的链接: v => e
        public List<int[]> v2e = new List<int[]>();

        // 离散边界点
        public List<PointF> verticeArray { set; get; }
        // 骨架线段
        public List<LineSegment2DF> skeletonArray { set; get; }
        // 区域邻接表
        public List<int[]> regionAdj;
        // Images
        public Image<Gray, Byte> binaryImage { set; get; }

        #region CONSTRUCTION

        public Gridding(List<PointF> vertices)
        {
            verticeArray = vertices;
            cellArray = new ArrayList();
        }
        public Gridding() { verticeArray = new List<PointF>(); cellArray = new ArrayList(); }
        #endregion

        #region AUXILIARY FUNCTION

        private static int[][] MOVE = new int[8][] { new int[2] { 0, -1 }, new int[2] { 1, -1 }, new int[2] { 1, 0 }, new int[2] { 1, 1 }, new int[2] { 0, 1 }, new int[2] { -1, 1 }, new int[2] { -1, 0 }, new int[2] { -1, -1 } };

        PointF GetPixelNeighbor(PointF a, int dir_8)
        {
            a.X += MOVE[dir_8][0];
            a.Y += MOVE[dir_8][1];
            return a;
        }

        PointF ProportionalPoint(PointF a, PointF b, float lambda)
        {
            return new PointF(lambda * a.X + (1 - lambda) * b.X, lambda * a.Y + (1 - lambda) * b.Y);
        }

        Point PointF2Point(PointF p)
        {
            return new Point((int)p.X, (int)p.Y);
        }





        double GetAngle(PointF p1, PointF p2, PointF p3)
        {
            double l1 = Math.Sqrt((p3.X - p2.X) * (p3.X - p2.X) + (p3.Y - p2.Y) * (p3.Y - p2.Y));
            double l2 = Math.Sqrt((p2.X - p1.X) * (p2.X - p1.X) + (p2.Y - p1.Y) * (p2.Y - p1.Y));
            return ((p3.X - p2.X) * (p2.X - p1.X) + (p3.Y - p2.Y) * (p2.Y - p1.Y)) / l1 / l2;
        }


        /// <summary>
        /// 返回三角形的种类
        /// </summary>
        /// <param name="t">三角形</param>
        /// <returns>三角形内部边数</returns>
        public int TriangleType(Triangle2DF t)
        {
            int cnt = 0;
            var points = t.GetVertices();
            for (int i = 0; i < 3; i++)
            {
                if (InnerLine(new LineSegment2DF(points[i], points[(i + 1) % 3])))
                    cnt++;
            }
            return cnt;
        }

        /// <summary>
        /// Returns if a triangle is an Inner Triangle
        /// </summary>
        /// <param name="t"></param>
        /// <returns></returns>
        public bool InnerTriangle(Triangle2DF t)
        {

            for (float l = 0.1F; l < 1; l += 0.1F)
                if (!InnerPoint(ProportionalPoint(t.V0, t.V1, l)))
                    return false;
            for (float l = 0.1F; l < 1; l += 0.1F)
                if (!InnerPoint(ProportionalPoint(t.V2, t.V1, l)))
                    return false;
            for (float l = 0.1F; l < 1; l += 0.1F)
                if (!InnerPoint(ProportionalPoint(t.V0, t.V2, l)))
                    return false;
            return true;
        }
        public bool InnerPoint(PointF p)
        {
            return (InImageRange(p, binaryImage) && binaryImage[PointF2Point(p)].Intensity > 0);
        }

        bool InImageRange(PointF p, IImage img)
        {
            return (p.X >= 0 && p.X <= img.Size.Width) && (p.Y >= 0 && p.Y <= img.Size.Height);
        }

        bool InnerLine(PointF p, PointF q, bool vertix = false)
        {
            if (Math.Min(Math.Abs(p.X - q.X), Math.Abs(p.Y - q.Y)) > 10)
            {
                float begin = vertix ? 0F : 0.3F;
                float end = 1 - begin;
                for (float k = begin; k < end; k += 0.3F)
                {
                    if (!InnerPoint(ProportionalPoint(p, q, k)))
                        return false;
                }
                return true;
            }
            else
            {
                return InnerPoint(new PointF((p.X + q.X) / 2, (p.Y + q.Y) / 2));
            }
        }



        bool InnerLine(LineSegment2DF l)
        {
            return InnerLine(l.P1, l.P2);
        }

        Image<Rgb, Byte> show = new Image<Rgb, byte>(800, 800);
        void VisualizeObject(IConvexPolygonF obj, ImageBox imageBox, Rgb rgb)
        {

            show.Draw(obj, rgb, 1);
            imageBox.Image = show;
        }
        void VisualizeObject(CircleF obj, ImageBox imageBox, Rgb rgb)
        {
            show.Draw(obj, rgb, 1);
            imageBox.Image = show;
        }

        void VisualizeObject(CircleF obj, string str, ImageBox imageBox, Rgb rgb)
        {
            MCvFont f = new MCvFont(Emgu.CV.CvEnum.FONT.CV_FONT_HERSHEY_COMPLEX, 0.4, 0.4); //Create the font
            show.Draw(str, ref f, PointF2Point(obj.Center), new Rgb(Color.Beige));
            show.Draw(obj, rgb, 1);
            imageBox.Image = show;
        }
        
        void VisualizeObject(LineSegment2DF obj, ImageBox imageBox, Rgb rgb)
        {
            show.Draw(obj, rgb, 1);
            imageBox.Image = show;
        }
        void VisualizeObject(IConvexPolygonF[] objs, ImageBox imageBox, Rgb rgb)
        {
            foreach (var obj in objs)
            {
                show.Draw(obj, rgb, 1);

            }
            imageBox.Image = show;
        }
        void VisualizeObject(CircleF[] objs, ImageBox imageBox, Rgb rgb)
        {
            foreach (var obj in objs)
            {
                show.Draw(obj, rgb, 1);

            }
            imageBox.Image = show;
        }
        void VisualizeObject(LineSegment2DF[] objs, ImageBox imageBox, Rgb rgb)
        {
            foreach (var obj in objs)
            {
                show.Draw(obj, rgb, 1);

            }
            imageBox.Image = show;
        }
        void VisualizeClear()
        {
            show = new Image<Rgb, byte>(800, 800);
        }
        #endregion

        #region GEOMETRICS

        /// <summary>
        /// if a triangle is clockwise
        /// </summary>
        /// <param name="t"></param>
        /// <returns></returns>
        public int ClockWise(Triangle2DF t)
        {
            Matrix<float> m1 = new Matrix<float>(3, 3);
            m1[0, 0] = t.V0.X;
            m1[0, 1] = t.V0.Y;
            m1[1, 0] = t.V1.X;
            m1[1, 1] = t.V1.Y;
            m1[2, 0] = t.V2.X;
            m1[2, 1] = t.V2.Y;
            m1[0, 2] = m1[1, 2] = m1[2, 2] = 1;
            if (m1.Det < -EPS)
                return -1;
            else if (m1.Det > EPS)
                return 1;
            else
                return 0;
        }

        private bool LineSegmentIntersect(MyPointF p1, MyPointF p2, MyPointF q1, MyPointF q2)
        {
            // vector dot mult
            int d1 = ClockWise(new Triangle2DF(q1, q2, p2)) * ClockWise(new Triangle2DF(q1, q2, p1));
            int d2 = ClockWise(new Triangle2DF(p1, p2, q1)) * ClockWise(new Triangle2DF(p1, p2, q2));
            return d1 < 0 && d2 < 0;
            //if (d1 < 0 && d2 < 0)
            //    return true;
            //else if (d1 > 0 || d2 > 0)
            //    return false;
            //else if (ClockWise(new Triangle2DF(q1, q2, p1)) != 0
            //    || ClockWise(new Triangle2DF(q1, q2, p2)) != 0)
            //    return false;
            //else
            //    return true;
        }

        public bool LineSegmentIntersect(PointF p1, PointF p2, PointF q1, PointF q2)
        {
            return LineSegmentIntersect(p1, p2, q1, q2);
        }

        public bool LineSegmentIntersect(LineSegment2DF l1, LineSegment2DF l2)
        {
            return LineSegmentIntersect(new MyPointF(l1.P1), new MyPointF(l1.P2), new MyPointF(l2.P1), new MyPointF(l2.P2));

        }

        public bool LineCoLine(LineSegment2DF l1, LineSegment2DF l2)
        {
            //           return (l1.P2.Y - l1.P1.Y) * (l2.P2.X - l2.P1.X) == (l2.P2.Y - l2.P1.Y) * (l1.P2.X - l1.P1.X);
            //return (Math.Abs(l1.Direction.Y - l2.Direction.Y) < EPS && Math.Abs(l1.Direction.X - l2.Direction.X) < EPS)
            //    || (Math.Abs(l1.Direction.Y + l2.Direction.Y) < EPS && Math.Abs(l1.Direction.X + l2.Direction.X) < EPS);
            return ClockWise(new Triangle2DF(l1.P1, l1.P2, l2.P1)) == 0 && ClockWise(new Triangle2DF(l1.P1, l1.P2, l2.P2)) == 0;
        }

        #endregion

        #region UTILITY

        public ProcedureState Compute()
        {

            // Exceptions
            if (verticeArray.Count == 0)
                return ProcedureState.FAIL;

            // Random Triangulation

            PointF[] curP = new PointF[3];
            int pointer = 0;
            curP[1] = (PointF)verticeArray[0];
            curP[2] = (PointF)verticeArray[1];
            for (int i = 2; i < verticeArray.Count; i++)
            {
                curP[pointer] = (PointF)verticeArray[i];
                cellArray.Add(new Triangle2DF(curP[0], curP[1], curP[2]));
                pointer = (pointer + 1) % 3;
            }
            return ProcedureState.SUCCESS;
        }

        /// <summary>
        /// Get Vertices From Canny Contour
        /// </summary>
        /// <param name="contour">Contour Image, 255true</param>
        /// <param name="scale">Discretelization Parameter, length of segments</param>
        /// <returns></returns>
        public List<PointF> GetVerticesFromContour(Image<Gray, Byte> contour, int scale = 10)
        {
            verticeArray = new List<PointF>();

            // procedure points
            PointF s = new PointF();
            PointF b, c;
            int i;
            // 从1开始，略去边上一圈，避免麻烦
            for (i = 1; i < contour.Height; i++)
            {
                bool BREAK = false;
                for (int j = 1; j < contour.Width; j++)
                {
                    if (contour[i, j].Intensity > 10)
                    {
                        s = new PointF(j, i);
                        BREAK = true;
                        break;
                    }
                }
                if (BREAK) break;
            }

            // c is set to the start point
            c = s;
            // b is the point right to the WEST of c
            b = c; b.X--;
            // step
            int count = 0;
            while (c != s || verticeArray.Count == 0)
            {
                if (count % scale == 0)
                    verticeArray.Add(c);
                count++;
                // repeatition once
                int ib;
                for (ib = 0; ib < 8; ib++)
                    if (GetPixelNeighbor(c, ib) == b) break;
                if (ib == 8) throw (new Exception("Contour Discretelization Error"));
                i = 0;
                for (i = 0; i < 8; i++)
                {
                    if (GetPixelNeighbor(c, (ib + i) % 8).ToString().Contains('-') || GetPixelNeighbor(c, (ib + i) % 8).X >= contour.Width || GetPixelNeighbor(c, (ib + i) % 8).Y >= contour.Height)
                        throw (new Exception("Neighbor Search Out_of_Range"));
                    if (contour[PointF2Point(GetPixelNeighbor(c, (ib + i) % 8))].Intensity > 10) break;
                }
                if (i == 8) throw (new Exception("Contour Discretelization Error"));
                c = GetPixelNeighbor(c, (ib + i) % 8);
                b = GetPixelNeighbor(c, (ib + i + 7) % 8);
            }

            return verticeArray;
        }

        /// <summary>
        /// check if the point is within the circumcircle of triangle
        /// </summary>
        /// <param name="p"></param>
        /// <param name="t"></param>
        /// <returns></returns>
        public bool PointInCircumcircle(PointF p, Triangle2DF t)
        {
            /*
             * 
             *  Time Reduction Needed:
             *  in fact there are two 3-degree determinant, so an explicit expression is availble
             * 
             */
            // regulate points of t in clock wise order
            Matrix<float> m1 = new Matrix<float>(3, 3);
            m1[0, 0] = t.V0.X;
            m1[0, 1] = t.V0.Y;
            m1[1, 0] = t.V1.X;
            m1[1, 1] = t.V1.Y;
            m1[2, 0] = t.V2.X;
            m1[2, 1] = t.V2.Y;
            m1[0, 2] = m1[1, 2] = m1[2, 2] = 1;
            if (m1.Det < 0)
            {
                PointF temp = t.V0;
                t.V0 = t.V1;
                t.V1 = temp;
            }
            // compute
            Matrix<float> m2 = new Matrix<float>(4, 4);
            m2[0, 0] = t.V0.X; m2[0, 1] = t.V0.Y; m2[0, 2] = (t.V0.X * t.V0.X + t.V0.Y * t.V0.Y); m2[0, 3] = 1;
            m2[1, 0] = t.V1.X; m2[1, 1] = t.V1.Y; m2[1, 2] = (t.V1.X * t.V1.X + t.V1.Y * t.V1.Y); m2[1, 3] = 1;
            m2[2, 0] = t.V2.X; m2[2, 1] = t.V2.Y; m2[2, 2] = (t.V2.X * t.V2.X + t.V2.Y * t.V2.Y); m2[2, 3] = 1;
            m2[3, 0] = p.X; m2[3, 1] = p.Y; m2[3, 2] = (p.X * p.X + p.Y * p.Y); m2[3, 3] = 1;

            return m2.Det >= 0;
        }

        public ArrayList GetGriddingFromVertices(ImageBox imgBox)
        {
            // get lines
            List<LineSegment2DF> lineArray = new List<LineSegment2DF>();
            for (int i = 0; i < verticeArray.Count; i++)
            {
                for (int j = i + 2; j < verticeArray.Count; j++)
                {
                    if ((j + 1) % verticeArray.Count == i) continue;
                    lineArray.Add(new LineSegment2DF(verticeArray[i], verticeArray[j]));
                }
            }
            IEnumerable<LineSegment2DF> sortAscending =
                from line in lineArray
                orderby line.Length ascending
                select line;

            // get line segments
            int lineCount = 0;
            lineSegmentArray = new List<LineSegment2DF>();
            // contour linesegments
                // point adj table
                bool[,] adj = new bool[verticeArray.Count, verticeArray.Count];
            for (int j = 0; j < verticeArray.Count - 1; j++)
            {
                lineSegmentArray.Add(new LineSegment2DF(verticeArray[j], verticeArray[j + 1]));
                adj[j, j + 1] = adj[j + 1, j] = true;
            }
            lineSegmentArray.Add(new LineSegment2DF(verticeArray[0], verticeArray[verticeArray.Count - 1]));
            adj[0, verticeArray.Count - 1] = adj[verticeArray.Count - 1, 0 ] = true;


            if (DEBUG)
            {

                foreach (var ver in verticeArray)
                {
                    VisualizeObject(new CircleF(ver, 2), imgBox, new Rgb(Color.White));
                }
            }
            // get legal inner line segments
            foreach (var line in sortAscending)
            {
                if (DEBUG)
                {
                    VisualizeObject(line, imgBox, new Rgb(Color.Red));
                }
                
                if (!InnerLine(line))
                {
                    continue;
                }

                if (DEBUG)
                {
                    VisualizeObject(line, imgBox, new Rgb(Color.Blue));
                    MessageBox.Show("inner:blue");
                }

                bool cross = false; int j = 0;
                for (j = 0; j < lineSegmentArray.Count; j++)
                {
                    var provLine = lineSegmentArray[j];
                    if (LineSegmentIntersect(line, provLine) || LineCoLine(line, provLine))
                    {
                        cross = true;
                        break;
                    }
                }
                if (cross) continue;
                if (DEBUG)
                {
                    VisualizeObject(line, imgBox, new Rgb(Color.Yellow));
                    MessageBox.Show("!cross:Yellow");
                }


                // complete gridding
                if (lineCount++ >= (verticeArray.Count - 3))
                    break;
                lineSegmentArray.Add(line);
                for (int k = 0; k < verticeArray.Count; k++)
                {
                    for (int l = 0; l < verticeArray.Count; l++)
                    {
                        if (verticeArray[k] == line.P1 && verticeArray[l] == line.P2)
                            adj[l, k] = adj[k,l] = true;
                    }
                }
            }


            bool[] used = new bool[verticeArray.Count];
            int count = 0;
            // find the first line
            for (int i = 0; i < verticeArray.Count; i++)
            {
                if (adj[0, i])
                {
                    used[0] = used[i] = true;
                    break;
                }
            }


            // viz
            LineSegment2DF[] lss = new LineSegment2DF[ 10000 ];
            int p = 0;
            for ( int i = 0 ;i < verticeArray.Count; i ++ )
                for (int j = 0; j < verticeArray.Count; j++)
                {
                    if (adj[i, j])
                    {
                        lss[p++] = new LineSegment2DF(verticeArray[i], verticeArray[j]);
                    }

                }
            if (DEBUG)
            {
                VisualizeObject(lss, imgBox, new Rgb(Color.White));
            }
            

            bool find = false;
            int triangleCnt = 0;
            ////////////////////////////////  TODO:        队列优化搜索过程  ///////////////////////////////
            while (triangleCnt++ < verticeArray.Count - 2)
            {
                for (int i = 0; i < verticeArray.Count; i++)
                {
                    if (find) break;
                    if (used[i])
                        for (int j = 0; j < verticeArray.Count; j++)
                        {
                            if (find) break;
                            if (used[j] && adj[i, j])
                            {
                                //MessageBox.Show("i:" + i.ToString() + "j:" + j.ToString());


                                for (int k = 0; k < verticeArray.Count; k++)
                                {
                                    if (adj[i, k] && adj[j, k] && !(used[i] && used[j] && used[k]))
                                    {
                                        cellArray.Add(new Triangle2DF(verticeArray[i], verticeArray[j], verticeArray[k]));
                                        used[i] = used[j] = used[k] = true;


                                        if (DEBUG)
                                        {
                                            VisualizeObject(new Triangle2DF(verticeArray[i], verticeArray[j], verticeArray[k]), imgBox, new Rgb(Color.Red));

                                            for (int m = 0; m < verticeArray.Count; m++)
                                            {
                                                if (used[m])
                                                    VisualizeObject(new CircleF(PointF2Point(verticeArray[m]), 3), m.ToString(), imgBox, new Rgb(Color.Blue));
                                            }
                                        }

                                        break;
                                    }
                                }
                            }
                        }
                }
            }

            // Final: Get Skeleton !
            skeletonArray = new List<LineSegment2DF>();
            for (int i = 0; i < cellArray.Count; i++)
            {
                for (int j = i+1; j < cellArray.Count; j++)
                {
                    if (TriangleAdj((Triangle2DF)cellArray[i], (Triangle2DF)cellArray[j]))
                        skeletonArray.Add(new LineSegment2DF(((Triangle2DF)cellArray[i]).Centeroid, ((Triangle2DF)cellArray[j]).Centeroid));
                }
            }

            return cellArray;

   //         foreach (var line in lineArray)
   //         {
   //             if (!InnerLine(line))
   //                 continue;
   //             bool cross = false;
   //             foreach (var provLine in lineSegmentArray)
   //             {
   //                 if (LineSegmentIntersect(line, provLine))
   //                 {
   //                     cross = true;
   //                     break;
   //                 }
   //             }
   //             if (!cross)
   //                 lineSegmentArray.Add(line);
   //         }

   //         return cellArray;

   //         //// random triangulation
   //         //PointF[] curP = new PointF[3];
   //         //int pointer = 0;
   //         //curP[1] = (PointF)vertices[0];
   //         //curP[2] = (PointF)vertices[1];
   //         //for (int i = 2; i < vertices.Count; i++)
   //         //{
   //         //    curP[pointer] = (PointF)vertices[i];
   //         //    cellArray.Add(new Triangle2DF(curP[0], curP[1], curP[2]));
   //         //    pointer = (pointer + 1) % 3;
   //         //}

   //         count = 0;

   //         // EAR TRIMMING
   //         List<PointF> que = verticeArray;
   //         while (que.Count > 3)
   //         {
   //             // most extrude angle's cos
   //             double earAngle = 2;
   //             int earInd = -1;
   //             for (int i = 0; i < que.Count; i++)
   //                 if (GetAngle(que[(i + que.Count - 1) % que.Count], que[i], que[(i + 1) % que.Count]) < earAngle
   //                     && InnerTriangle(new Triangle2DF(que[(i + que.Count - 1) % que.Count], que[i], que[(i + 1) % que.Count]))
   //)
   //                 {
   //                     earAngle = GetAngle(que[(i + que.Count - 1) % que.Count], que[i], que[(i + 1) % que.Count]);
   //                     earInd = i;
   //                 }
   //             if (earInd < 0)
   //                 break;
   //             cellArray.Add(new Triangle2DF(que[(earInd + que.Count - 1) % que.Count], que[earInd], que[(earInd + 1) % que.Count]));
   //             que.RemoveAt((earInd + 1) % que.Count);
   //         }
   //         cellArray.Add(new Triangle2DF(que[0], que[1], que[2]));

   //         return cellArray;

   //         for (int i = 0; i < cellArray.Count; i++)
   //             RecursiveFlip(i);
   //         // // n^2 version
   //         //List<KeyValuePair<int, int>> que = new List<KeyValuePair<int,int>>();
   //         //for (int i = 0; i < cellArray.Count; i++)
   //         //    for (int j = i + 1; j < cellArray.Count; j++)
   //         //        if (TriangleAdj(((Triangle2DF)cellArray[i]), (Triangle2DF)cellArray[j]))
   //         //            que.Add(new KeyValuePair<int, int>(i, j));
   //         //List<int> neighbors = new List<int>();
   //         //while (que.Count > 0)
   //         //{
   //         //    // 取下标
   //         //    int n1 = que[0].Key; int n2 = que[0].Value;
   //         //    // 取三角形
   //         //    Triangle2DF t1 = (Triangle2DF)cellArray[n1];
   //         //    Triangle2DF t2 = (Triangle2DF)cellArray[n2];
   //         //    if (Flip(ref t1, ref t2))
   //         //    {
   //         //        // Get All Neighbors
   //         //        foreach (var link in que)
   //         //        {
   //         //            if (link.Key == n1 || link.Key == n2)
   //         //                neighbors.Add(link.Value);
   //         //            if (link.Value == n1 || link.Value == n2)
   //         //                neighbors.Add(link.Key);
   //         //        }
   //         //        for ( int i = 0; i < que.Count; i ++ ){
   //         //            KeyValuePair<int,int> link = que[i];
   //         //            if (link.Key == n1 || link.Key == n2)
   //         //            {
   //         //                if (TriangleAdj((Triangle2DF)cellArray[link.Value], (Triangle2DF)cellArray[n1]))
   //         //                    link.Value = n1;
   //         //                else
   //         //                    link.Value = n2;
   //         //            }
   //         //        }


   //         //    }
   //         //}

   //         return this.cellArray;
        }

        /// <summary>
        /// Try flip among two triangles. Return if Flip is performed
        /// </summary>
        /// <param name="a">reference of triangle a</param>
        /// <param name="b">reference of triangle b</param>
        /// <returns>If flipped</returns>
        bool Flip(ref Triangle2DF a, ref Triangle2DF b)
        {
            int i, k, j;
            // find single points in each of the triangle
            for (i = 0; i < 3; i++)
            {
                for (k = 0; k < 3; k++)
                    if (a.GetVertices()[i] == b.GetVertices()[k])
                        break;
                if (k == 3) break;
            }
            for (j = 0; j < 3; j++)
            {
                for (k = 0; k < 3; k++)
                    if (b.GetVertices()[j] == a.GetVertices()[k])
                        break;
                if (k == 3) break;
            }
            if (PointInCircumcircle(b.GetVertices()[j], a))
            {
                PointF pa = a.GetVertices()[i];
                PointF pb = b.GetVertices()[j];
                PointF p1 = a.GetVertices()[(i + 1) % 3];
                PointF p2 = a.GetVertices()[(i + 2) % 3];
                a = new Triangle2DF(p1, pa, p2);
                b = new Triangle2DF(p1, pb, p2);
                return true;
            }
            else
                return false;
        }

        /// <summary>
        /// 递归调用Flip
        /// </summary>
        /// <param name="i">三角形在cellArray中的Index</param>
        /// <returns>是否进行翻转</returns>
        bool RecursiveFlip(int ind)
        {
            bool change = false;
            for (int i = 0; i < cellArray.Count; i++)
            {
                Triangle2DF t1 = (Triangle2DF)cellArray[i];
                Triangle2DF t2 = (Triangle2DF)cellArray[ind];
                if (TriangleAdj(t1, t2) && Flip(ref t1, ref t2))
                {
                    cellArray[i] = t1;
                    cellArray[ind] = t2;
                    for (int j = 0; j < cellArray.Count; j++)
                        if (TriangleAdj((Triangle2DF)cellArray[j], t1))
                            RecursiveFlip(i);
                    change = true;
                }
            }
            return change;
        }

        /// <summary>
        /// check if two triangles are edge-Adjacent
        /// </summary>
        /// <param name="t1"></param>
        /// <param name="t2"></param>
        /// <returns></returns>
        bool TriangleAdj(Triangle2DF t1, Triangle2DF t2)
        {
            int i = 0;
            foreach (var p in t1.GetVertices())
                foreach (var pp in t2.GetVertices())
                {
                    if (p == pp) i++;
                }
            return i == 2;
        }




        #endregion

        #region AUXILIARY VISUALIZATION

        public Image<Rgb, Byte> ShowGridding()
        {
            Image<Rgb, Byte> canvas = new Image<Rgb, byte>(200, 200);

            foreach (PointF p in verticeArray)
            {
                canvas.Draw(new CircleF(p, 3), new Rgb(Color.Red), 1);
            }
            foreach (Triangle2DF t in cellArray)
            {
                canvas.Draw(t, new Rgb(Color.Blue), 1);
            }
            return canvas;
        }

        #endregion

    }
    class MyOCR
    {
        Gridding n = new Gridding();
    }
}
